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ABSTRACT 


This  thesis  addresses  the  problem  of  integrated  design  of  the  aircraft  plant  pa¬ 
rameters  and  of  the  corresponding  feedback  controller.  The  plant  parameters  are 
typically  the  sizes  of  the  control  surfaces  or  other  aerodynamical  surfaces  of  the 
aircraft.  The  approach  is  to  rewrite  the  aircraft  dynamic  requirements  as  linear 
matrix  inequalities  (LMI’s)  and  to  optimize  a  linear  cost  function  associated  with 
aircraft  plant  parameters,  while  meeting  the  LMI  constraints.  An  algorithm  using 
Matlab  and  LMI-Lab  has  been  developed.  This  algorithm  has  been  used  for  inte¬ 
grated  plant /controller  design  of  an  F-4  aircraft  at  five  different  flight  conditions. 
The  result  consists  of  a  set  of  controllers,  one  for  each  flight  condition,  and  a  single 
solution  for  the  optimal  sizes  of  the  F-4  stabilator  and  spoiler  control  surfaces. 
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I.  INTRODUCTION 


Control  systems  of  modern  aircraft  have  become  increasingly  complex.  That 
is  especially  true  in  the  case  of  fighter  aircraft  with  relaxed  stability  and  multiple 
control  surfaces  working  in  each  axis. 

The  control  system,  control  surfaces  included,  must  satisfy  all  the  performance 
requirements  (stability,  disturbance  rejection,  maneuverability  etc.),  while  not  causing 
excessive  drag,  weight  or  signatures.  It  is  thus  interesting  to  study  the  problem  of 
finding  the  smallest  size  of  some  parameters  describing  the  physical  configuration  of 
the  aircraft,  while  still  being  able  to  find  a  controller  that  satisfies  the  performance 
requirements. 

While  the  importance  of  integrated  plant  and  controller  design  has  been  recog¬ 
nized,  little  had  been  done  in  the  area  prior  to  [Ref.  1].  In  this  thesis  the  ideas  of 
integrating  a  minimization  of  the  control  surfaces  with  the  design  of  a  controller,  that 
were  developed  in  [Ref.  1],  are  developed  further. 

The  design  procedure  is  based  on  a  set  of  requirements  that  are  to  be  fulfilled  at 
a  number  of  different  flight  conditions  and  a  cost  function  for  the  plant  parameters, 
e.g.  the  size  of  the  control  surfaces.  The  result  consists  of  a  set  of  controllers,  one 
for  each  flight  condition,  and  one  single  optimal  solution  for  the  size  of  the  plant 
parameters.  The  proposed  design  procedure  is  still  in  an  early  prototype  phase,  but 
shows  promise  of  becoming  a  useful  tool. 

The  thesis  is  organized  as  follows.  Chapter  II  contains  the  formulation  of  the 
treated  problem  as  an  Hoo  control  problem,  and  the  parametrization  of  the  plant.  In 
Chapter  III  the  algorithm  from  Chapter  II  is  applied  to  the  longitudinal  dynamics 
of  an  F-4  fighter  aircraft.  The  plant  parameters  are  chosen  to  be  the  sizes  of  the 
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moving  horizontal  tail  and  of  the  spoiler.  Chapter  II  also  contains  the  results  from 
the  plant  optimization  and  the  associated  controllers.  The  program  code  in  Matlab 
and  LMI-Lab  that  is  used  to  solve  this  problem  is  included  in  the  appendix. 
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II.  PROBLEM  STATEMENT 


The  task  of  designing  the  control  system  of  an  aircraft  has  become  increasingly 
complex  in  recent  time.  That  is  especially  true  in  the  case  of  modern  fighter  aircraft 
with  relaxed  stability  and  many  control  surfaces  working  in  each  axis.  The  control 
system  must,  of  course,  meet  all  the  imposed  performance  requirements,  but  as  ex¬ 
cessive  control  capability  costs  weight,  space,  money,  signature  etc.,  it  is  desireable 
to  optimize  the  control  system. 

The  problem  to  be  addressed  is  to  find  the  minimal  size  of  some  plant  parameters 
(typically  the  size  of  the  control  surfaces  or  other  aerodynamical  surfaces),  while 
finding  a  feedback  controller  that  meets  the  stated  control  requirements.  The  design 
requirements  are  given  for  several  flight  conditions  and  the  designed  combination  of 
plant  parameter  and  controller  must,  of  course,  work  for  all  of  the  flight  conditions. 
The  result  of  the  optimization  is  one  set  of  values  for  the  plant  parameters  and  a 
set  of  controllers,  one  for  each  flight  condition.  The  objective  is  to  use  this  set  of 
controllers  to  design  a  gain  scheduled  controller,  that  works  for  all  the  given  flight 
conditions  (and  hopefully  in  between  those  given  conditions). 

When  the  sizes  of  the  surfaces  are  varied,  the  dynamics  of  the  aircraft  changes. 
It  is  therefore  natural  to  integrate  the  minimization  of  a  cost  function  associated  with 
the  plant  parameters  with  the  controller  design.  The  controller  is  designed  using  H^o 
methods,  which  are  summarized  in  Section  A.  The  problem  of  minimizing  the  plant 
parameters,  while  still  being  able  to  find  a  feasible  controller,  is  treated  in  Section  B. 
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A.  Hoo  CONTROL 


Figure  2.1:  Feedback  interconnection  of  plant  Q  and  controller  fC. 

A  general  feedback  system  with  exogenous  input  w  and  output  to  be  controlled 
z  is  shown  in  Figure  2.1.  In  this  context  we  will  only  consider  linear,  time  invariant 
systems.  Here  u;  is  a  reference  signal  that  includes  disturbance  inputs.  We  wish 
the  system  to  reject  the  disturbance  and  track  the  reference  signal,  i.e.  we  want 
a  controller  that  makes  the  dilference  between  the  reference  input  and  the  actual 
output  &s  small  as  possible,  while  not  expending  too  much  control  energy.  This 
is  accomplished  by  making  z,  which  is  a  combination  of  tracking  errors  and  the 
control  inputs,  as  small  as  possible.  Let  Tzw{G,  K)  denote  the  closed  loop  transfer 
fuction  from  the  input  w  to  the  output  z.  The  solution  to  the  Hoo  control  problem 
is  the  controller  K  that  minimizes  the  infinity  norm  of  the  transfer  function,  i.e. 
\\TzUG,K)\\oo. 

The  infinity  norm  is  a  generalization  of  the  maximal  gain  and  is  defined  as  the 
supremum  over  all  frequencies  of  the  maximal  singular  value  of  the  transfer  function, 
i.e. 

=  sup{a(r^„,(ia;))} 

LU 

Although  the  above  definition  is  given  in  terms  of  the  frequency  domain,  state- 
space  methods  are  very  useful  in  Hoo  synthesis.  A  state-space  realization  of  the  plant 
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depicted  in  Figure  2.1  with  state  feedback  can  be  written  as 


Ax  -t-  B\w  -f  B2U 

(2.1) 

Cx  Du 

(2.2) 

X 

(2.3) 

Where  x  €  jR",  ^  u  G  iZ™  and  w  €  iZ®. 

To  ensure  the  existence  of  a  solution  to  the  /Zoo  control  problem,  we  have  to 
assume  that  {C2,  A,  B2)  is  stabilizable  and  detectable,  and  that  D  has  linearly  inde¬ 
pendant  columns.  The  inputs  and  outputs,  i.e.  w  and  z,  are  usually  scaled  in  such  a 
way  that  the  design  requirements  are  satisfied  if  /^)||oo  <  1- 

It  can  be  shown  [Ref.  2]  that  there  exists  a  controller  K  such  that  \\Tzw{G,  /f  )|loo  < 
1  if  and  only  if  there  exist  W  and  Y,  with  Y  symmetric  and  positive  definite,  such 
that 

'  AY +  YA^  +  B2W +  W^BJ  Bi  (CY  +  DW)'^' 

R{W,Y)=  Bj  -1  0  <0  (2.4) 

CY^DW  0  -I 

If  the  linear  (or  rather  affine)  matrix  inequality  (2.4)  has  a  solution,  the  con¬ 
troller  K  =  WY~^  solves  the  //oo  control  problem. 

B.  PLANT  PARAMETRIZATION 

We  consider  a  dynamical  system  at  several  given  operating  conditions,  or  more 
specifically  an  aircraft  with  performance  requirements  given  for  several  flight  condi¬ 
tions.  At  each  of  the  flight  conditions  we  use  a  linear  model  to  describe  the  dynamics 
of  the  aircraft.  The  linear  models  are  obtained  by  linearizing  the  general  non-linear 
six  degree  of  freedom  equations  of  motion  for  each  flight  condition.  The  linearized 
models  are  different  for  each  flight  condition  and  with  i,  i  =  1,2,  ...Z,  denoting  the 
respective  flight  conditions  we  have 
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(2.5) 

(2.6) 


X  =  A'x  +  Blw  +  B2U 
z  =  C‘x  +  D^u 


Where  x  €  i?”  is  the  state-vector  representing  perturbations  around  the  trimmed 
condition,  w  €  is  the  external  input,  u  €  R”^  is  the  control  input  and  z  €  R’’  is 


the  ouput  to  be  regulated. 

We  now  parametrize  the  plant  with  the  plant  parameter  (  —  [Ci? Cz? •••> 


where  C  is  normalized  such  that  0  <  Cj  <  1-  Associated  with  C  is- a  linear  cost  function 
J  =  c^C,  which  we  wish  to  minimize.  The  and  B2  matrices,  describing  the 

plant  are  parametrized  as  the  first  terms  of  the  Taylor  expansion  around  the  point 

^  =  0  as 

-4’(c)  = 

j  =  l 

(2.7) 

Bid)  =  BU  +  'Z.Bidi 

(2.8) 

Bid)  =  + 

(2.9) 

where 

A’  -  [^'^'1 
[^0Jc=o 

(2.10) 

D*  _  A 

“  dC 

L  J  c=o 

(2.11) 

n-  - 

~  dC 

L  J  c=o 

(2.12) 

The  C’  and  £)'  matrices  could  also  depend  on  but  that  case  is 

not  considered 

here. 

For  each  of  the  flight  conditions  we  can  find  a  steady  state  Hoo  controller  by 

solving  (2.4)  for  W'  and  Y\  However,  the  plant  matrices  now  depend  on  (,  hence 
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■  A*(C)y  +  +  B\{C,)W^  +  (C)  Bi(C)  (C*r‘  +  D^W^f  ' 

R(W‘,Y\0=  5^(0  -I  0  <0 

c^'y’  +  0  -I 

(2.13) 

Our  problem  is  thus 
min  J  = 

s.t.  (y‘,iy’,c)  €  $ 
where  the  feasible  set  $  is  defined  as 

$  =  {w^  e  F’  G  c :  F‘  >  0, y  =  y^, 0  <  0  <  i, R(vy',  y , c)  <  o} 

C.  ALGORITHM 

Following  the  results  in  [Ref.  1]  the  optimization  problem  in  W*,  Y'  and  ^  is 
solved  as  follows.  We  observe  that  R,  for  a  fixed  is  affine  in  the  two  remaining 
variables  W'  and  Y\  We  further  observe  that  R,  for  fixed  F*  and  W',  is  affine  in  the 
sole  remaining  variable 

Hence  the  approach  is  to: 

1.  Fix  C  —  Co  and  find  W'  and  F’  >  0,  i  =  1,2, ...,/  such  that  R(iy  ,F‘,Co)  <  0 

2.  With  VF*  and  F'  from  step  1,  find  the  C  with  Cj  >  0  that  minimizes  J  =  c^C 
s.t.  R(W',F‘,C)  <0 

3.  Go  to  step  1.  until  exit  criterion  is  satisfied. 

The  problem  to  be  solved  in  step  1  is  to  solve  I  LMI’s  for  the  matrices  W' 
and  F‘.  The  /  LMFs  are,  however,  independent  of  each  other  and  can  be  solved  in 
sequence.  The  problem  in  step  2  is  that  of  solving  one  single  LMI  for  the  minimal 
plant  parameter  C-  Here  all  flight  conditions  must  be  treated  simultaneously,  as  the 
same  physical  aircraft  must  be  used  at  all  the  flight  conditions. 
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The  solution  of  LMI’s  is  supported  by  software  packages  such  as  LMI-Lab  [Ref. 
6],  which  is  used  in  the  example  given  in  Chapter  III.  While  the  approach  to  formu¬ 
late  the  problem  as  two  LMI’s  solved  in  sequence  makes  the  calculations  manageble, 
there  is  unfortunately  no  guarantee  that  the  resulting  algorithm  will  converge  to  a 
stationary  point,  [Ref.  3].  Neither  is  it  clear  whether  or  not  the  feasible  set  $  is 
convex.  Convexity  is,  of  course,  of  utmost  importance  in  every  optimization  problem, 
as  when  a  convex  function  is  minimized  over  a  convex  set,  all  minimas  are  global. 

Chapter  III  describes  the  application  of  this  algorithm  to  the  longitudinal  dy¬ 
namics  of  an  F-4  lighter  aircraft. 
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III.  LONGITUDINAL  CONTROL  OF  F-4 


The  specific  example  concerns  the  control  of  the  longitudinal  dynamics  of  an 
F-4C  fighter  aircraft.  The  objective  is  to  size  the  two  control  surfaces,  elevator  and 
spoiler,  and  to  design  a  controller  for  elevator,  spoiler  and  thrust  for  each  of  the 
different  flight  conditions.  Five  flight  conditions,  three  subsonic  and  two  supersonic, 
with  data  from  [Ref.  4]  are  used;  see  tables  3.1  and  3.2. 

The  problem  formulation  is  basically  the  same  as  in  [Ref.  1],  i.e.  to  minimize 
the  sizes  of  the  elevator  and  the  spoiler,  while  finding  a  controller  that  satisfies  a 
set  of  performance  requirements  for  each  flight  condition.  The  main  difference  is 
that  thrust  control  is  included  and  that  multiple  flight  conditions  are  treated.  The 
numerical  algorithm,  which  uses  LMI-Lab,  is  completely  different  from  the  one  used 
in  [Ref.  1]. 

In  section  A  the  equations  of  motion  are  parametrized  according  to  the  chosen 
parameters,  elevator  and  spoiler.  Section  B  treats  the  performance  requirements  and 
how  they  are  expressed  in  terms  of  /foo  control.  Sections  C-E  contain  the  results  from 
the  performed  calculations. 

Nomenclature: 

Xa  Aerodynamic  derivative  with  respect  to  states  or  controls  in  X-direction 

Za  Aerodynamic  derivative  with  respect  to  states  or  controls  in  Z-direction 

Ma  Aerodynamic  derivative  with  respect  to  states  or  controls  in  M-direction 

ly  Moment  of  inertia  around  the  y-axis 
m  Mass 

U  Trimmed  velocity  in  x-direction 
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Qo  Trimmed  angle  of  attack 

I  Distance  from  center  of  gravity  to  elevator 

6e  Elevator  deflection 

Ss  Spoiler  deflection 

6t  Thrust  perturbation  from  equilibrium 

(7a  Disturbance  in  angle  of  attack,  a 

<7^  Disturbance  in  flight  path  angle,  7 

A.  EQUATIONS  OF  MOTION  AND  AERODYNAMICAL  DATA 

The  equations  describing  the  dynamics  of  an  aircraft  in  six  degrees  of  freedom 
consist  of  a  set  of  coupled  non-linear  differential  equations.  These  equations  are 
conveniently  formulated  in  a  body  fixed  coordinate  system.  In  the  design  of  the 
control  system,  we  generally  consider  the  set  of  linearized  equations  describing  the 
dynamics  around  a  certain  nominal  flight  condition. 

In  this  section  we  use  a  model  with  four  linear  differential  equations  to  describe 
the  longitudinal  dynamics.  In  general  the  longitudinal  dynamics  of  an  aircraft  can  be 
divided  into  two  modes:  the  short  period  and  the  long  period.  From  [Ref.  4]  we  obtain 
aerodynamical  data  for  an  F-4  fighter  aircraft  operating  at  five  flight  conditions:  three 
subsonic  and  two  supersonic. 

The  plant  parameters  to  be  optimized  are  the  size  of  the  moving  horizontal 
tail  and  the  size  of  a  spoiler.  These  sizes  are  represented  by  the  two  coefficients 
Zs  and  Zs,  which  give  the  vertical  acceleration  per  radian  of  elevator  and  spoiler 
deflection  respectively.  Hence,  the  plant  parameter  is  C  = 

aerodynamical  data  is  adapted  to  fit  that  parametrization,  by  first  calculating  the 
aerodynamical  derivatives  for  the  condition  C  ~  ^  then  calculating  all  the  matrices 
needed  in  the  equations  of  motion. 
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TABLE  3.1:  Aircraft  Data  (Ref.  4) 


Flight  Condition 

(2) 

■M 

(5) 

Altitude  {ft) 

15,000 

35,000 

35,000 

35,000 

45,000 

Mach 

0.6 

0.9 

1.2 

1.5 

V{ft/s) 

952 

584 

876 

1167 

1452 

126 

489 

1210 

1210 

1210 

1210 

1210 

Ix{slug  •  fP) 

24970 

27360 

25040 

24970 

25040 

Iy(slug  ■  fp) 

122190 

122190 

122190 

122190 

122190 

Iz{slug  •  fP) 

13970 

139800 

IB 

139730 

hz{slug  •  fP) 

1175 

-16432 

-3033 

-1030 

-3033 

9.4 

1.6 

2.6 

The  spoiler  is  an  add-on  to  the  data  from  [Ref.  4]  and  is  modeled  in  the  same 
manner  as  in  [Ref.  1],  i.e.  as  giving  a  maximal  acceleration  of  0.2g  in  the  2-direction 
and  a  corresponding  acceleration  in  the  x-direction  of  —0.02g. 

We  need  to  find  out  how  the  varying  of  the  plant  parameters  affect  the  aerody¬ 
namical  derivatives  describing  the  aircraft  dynamics.  Changing  the  size  of  the  moving 
tail  affects  the  dynamics  of  the  aircraft  drastically.  As  the  tail  is  reduced  the  stability 
as  well  as  the  damping  is  reduced.  The  other  plant  parameter,  the  spoiler,  does, 
however,  not  affect  any  coefficients  other  than  and  Zg.  The  division  of  the  aero¬ 
dynamic  derivatives  into  one  part  that  depends  on  the  wing  and  body  of  the  aircrft 
and  one  part  that  depends  on  the  tail  is  done  according  to  [Ref.  5]  as 


ZqixjjI)  -  Zot  Z^ 

Zqwb  ~ 

I 

^awb 

^awb  ~  Miy  ~  Zg 

7/ 


■^qwb  —  ^^q 


Pm 

1^ 


Zs 
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(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 


TABLE  3.2:  Nominal  Dimensional  Derivatives  (Ref.  4) 


Flight  Condition 

(1) 

(2) 

(3) 

(4) 

(5) 

Xu{s-^) 

-0.0215 

-0.0176 

-0.0123 

-0.0136 

-0.0072 

-5.93 

-24.25 

-7.43 

-16.53 

-29.98 

Zu{s-^) 

-1.170 

-0.337 

-0.571 

-0.747 

-0.515 

Zaift/S^) 

-1103.2 

-162.2 

-475.4 

-848.3 

-716.7 

Zaift/s) 

-2.00 

-0.591 

-1.01 

-1.24 

-0.519 

-6.00 

-1.82 

-2.89 

-4.09 

-2.24 

Muift 

-0.0044 

-0.0001 

-0.0028 

0.0022 

0.0025 

Ma{s-'^) 

-17.00 

-1.927 

-7.887 

-29.03 

-28.94 

Mais~l) 

-0.457 

-0.141 

-0.234 

-0.288 

-0.122 

-0.993 

-0.307 

-0.487 

-0.746 

-0.488 

Xs(Stls‘) 

0.00 

-0.01 

0.00 

-0.01 

0.00 

-107.0 

-20.98 

-49.65 

-90.44 

-70.67 

Ms{s-^) 

-25.00 

-4.90 

-11.40 

-20.70 

-16.00 

-0.919 

-0.919 

-0.919 

-0.919 

-0.919 

z.ifth’) 

9.19 

9.19 

9.19 

9.19 

9.19 

=  Ma-jjjZs  (3.6) 

The  resulting  wing/body  aerodynamical  derivatives  are  shown  in  Table  3.3.  A 
comparison  with  the  values  in  Table  3.2  shows  that  the  removal  of  the  horizontal  tail 
causes  the  static  stability  to  decrease.  In  the  three  subsonic  conditions  the  aircraft 
is  now  unstable  {Mo,  >  0).  The  damping  has  also  decreased,  i.e.  Ma  +  Mg  is  now 
closer  to  zero.  The  fact  that  Ma  is  positive  is  unusual,  but  as  it  is  usually  difficult 
to  separate  Ma  and  Mg,  and  they  still  add  up  to  a  negative  number,  there  is  no  need 
for  concern. 

For  the  longitudinal  dynamics  we  use  the  linearized  equations  of  motion  [Ref.  4] 
with  the  four  state  variables  given  as  [u/U,a,q,d]^ .  These  states  represent  perturba¬ 
tions  around  the  trimmed  values  [U,  oq,  0,  By  using  the  dimensional  aerodynamic 
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TABLE  3.3:  Wing/Body  Dimensional  Derivatives 


Flight  Condition 

(1) 

(2) 

(3) 

(4) 

(5) 

^«(s-*) 

-0.0215 

-0.0176 

-0.0123 

-0.0136 

-0.0072 

XaUth'') 

-5.93 

-24.25 

-7.43 

-16.53 

-29.98 

Zu{s-^) 

-1.170 

-0.337 

-0.571 

-0.747 

-0.515 

Zaiftls^) 

-996.2 

-141.2 

-425.75 

-757.9 

-646.0 

Zaifijs) 

0.652 

0.256 

0.318 

0.551 

0.594 

-3.35 

-0.973 

-1.56 

-2.30 

-1.13 

-0.0044 

-0.0001 

-0.0028 

0.0022 

0.0025 

8.00 

2.97 

3.52 

-8.33 

-12.94 

0.163 

0.057 

0.071 

0.122 

0.130 

-0.373 

-0.109 

-0.182 

-0.336 

-0.236 

derivatives  the  equations  of  motion  for  the  longitudinal  case  assume  the  form 


rj  =  +  B{^w  +  B. 


2n 


U 


(3.7) 


where  w  =  [ore/?  7re/]^  is  the  external  input  and  u  =  is  the  commanded 

values  of  elevator  deflection,  spoiler  deflection  and  thrust. 

With  aerodynamical  data  for  the  respective  flight  condition,  the  plant  matrices 
are  given  as 


■  u 

0 

0 

0  ■ 

0  u 
0  - 

-z^  0 

-Ma  1 

0 

0 

.  0 

0 

0 

1 . 

■  UX, 

,  X, 

0 

—gcosOo 

UZu 

Za 

u  +  z, 

—g  sin  ^0 

UM^ 

M, 

0 

0 

0 

0 

1 

■  0  ■ 
Za  0 
Ma  0 
0  0 


(3.8) 


(3.9) 


(3.10) 
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0  X  1 

Zg  0 
0  0 
0  0  0 


(3.11) 


(3.12) 


(3.13) 


U  =  Sg 

Ji  . 

We  now  write  the  A*,  B\^  and  matrices  as 

2 


(3.14) 


.4‘(C) 

=  .4‘  +  2/l'C 

J  =  1 

(3.15) 

B{{0 

i=i 

(3.16) 

B\{Q 

=  fljo+EfiiA 

(3.17) 

This  is  the  equivalent  of  taking  the  first  two  terms  in  the  Taylor  series  expansion 
around  the  point  (  =  [0,0]^.  Using  the  matrices  defined  in  (3.8-3.11)  we  have 


dCi  dZs  ^  dZs  ”  ^  dZs  ’ 

dC2  dZg  ™  '  dZg  "  ”  dZg  ’ 


dB\a  _  dB 
5Ci  “  dZ 


^7 

^Snom  V  Q2^  -^In  '  v-^n/  j^Snom 


(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 
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dB\Q 

9K 

d(^2 

dZ, 

[liv 

dBi, 

dBio 

dCi 

dZs 

dBio 

_  ^^20 

r7 

^  2"  ^snom 


+/rr’^>7 

'  dZ.  ''”'  dZ.  '^' 


V  ^2n  ’  \^n)  dZg 

~  V  ^2n  f  \^n)  Q2i 


(3.23) 

(3.24) 

(3.25) 

(3.26) 


Where  all  the  matrices  are  evaluated  at  [Zg,  Z*]  =  [0,0].  We  now  calculate  the 
needed  derivatives. 


dZs 


■  0  0  0  0 

10  1  0  0 

u{u-Zg,y  0  M<i  +  ^(c/-z^)  0  0 

0  “o  0  0 


(3.27) 


dZs 


(3.28) 


0  0  0  0 

0  1  1/t/  0 

n  hn  l^m  n 

^  ly  lyU  ^ 

0  0  0  0 


(3.29) 


^=0 

dZ, 


(3.30) 


^  0 
0  0 


(3.31) 


(3.32) 


0  0  0 

1  0  0 

^00 

h 

0  0  0 


(3.33) 
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(3.34) 


■  0  xjz,  0  ■ 

0  1  0 
dZs  0  0  0 

.0  0  0  . 

B.  DESIGN  OBJECTIVES  AND  AUGMENTED  SYSTEM 

The  design  objectives  follow  those  outlined  in  [Ref.  1]: 

•  Step  response  -  There  should  be  no  steady  state  error  in  the  response  to  a  step 
input  in  angle  of  attack,  a,  or  flight  path  angle,  'y  =  9  —  a. 

•  Closed  loop  must  be  internally  stable. 

•  Gust  disturbaces  with  a  magnitude  of  5  ft/s  in  the  z-direction  should  not  yield 
elevator  deflection  greater  than  20  degrees,  spoiler  deflection  greater  than  40 
degrees  or  errors  in  angle  of  attack  greater  than  1.5  degrees. 

•  A  flight  path  angle  of  1  degree  must  be  generated  using  elevator  and  spoiler 
deflection  of  less  than  20  and  40  degrees,  respectively. 


To  accommodate  the  steady  state  tracking  requirement,  the  plant  is  augmented 
with  the  integrals  of  a  and  7.  The  augmented  system  has  six  states,  x  =  [u/t/,  a,  q,  0,  a/s,  7/s]^. 
The  parametrized  augmented  system  is  thus 


j=l  j=l  j+1 

=  Cixa  +  Diu  (3.36) 


with  the  respective  matrices  given  as 


0  0  ■ 

Ai  0  0 

0  0 

0  0 

0  0  0  0 

0  10  0. 


(3.37) 
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(3.38) 

(3.39) 

(3.40) 

(3.41) 

(3.42) 

(3.43) 

(3.44) 

(3.45) 


(3.46) 


•®22 


0 

0 


0 

0 


C.  CALCULATIONS 

The  calculations  were  done  in  Matlab.  The  code  is  listed  in  the  appendix.  The 
linear  matrix  inequalities  are  solved  by  using  LMI-Lab,  [Ref.  6],  which  is  a  recently 
developed  toolbox  to  Matlab.  W‘  and  i.e.  the  controllers,  are  obtained  from  the 
function  ’’feasp”  and  the  minimization  of  (  is  done  by  using  ’’linobj”.  The  function 
’’feasp”  sometimes,  however,  returns  a  solution  that  is  not  at  all  negative  definite. 
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That  problem  is  solved  by  testing  the  calculated  solution  and,  when  it  is  infeasible, 
discarding  it  in  favor  of  the  solution  from  the  previous  iteration.  When  ’’feasp”  is 
again  called  in  the  next  iteration,  it  usually  gives  a  feasible  solution. 

The  code  is  fairly  fast  -  a  run  takes  about  10  minutes  on  a  Sparc  10. 

The  calculations  were  performed  with  the  following  values: 

■  1  ■ 

"=[l_ 

<T.y  =  1  degree 
a  a  =  5/Ui 
Semax  =  20  degrees 
^smax  =  40  degrees 
^tmax  —  0.05^ 

D.  OPTIMIZATION  RESULTS 

As  we  are  uncertain  about  whether  or  not  the  feasible  set  $  is  convex,  the  con¬ 
vergence  properties  of  the  algorithm  were  tested  to  see  if  the  solution  would  converge 
to  a  single  optimal  point  or  not. 

The  convergence  properties  were  tested  by  running  the  algorithm  with  different 
initial  conditions.  First  we  varied  the  initial  value  of  the  elevator,  i.e.  (unitiaii  between 
0.1  and  1.5,  while  holding  the  initial  value  of  the  spoiler  as  the  nominal,  i.e.  C2imtiai  = 
1.  In  3.1  we  see  the  optimization  history  for  each  iteration  of  the  two  plant  parameters. 
The  solution  tends  to  converge  to  smaller  values,  the  smaller  the  initial  value  is.  The 
same  results  are  shown  in  3.2,  where  the  two  values  i^i  and  (2  are  plotted  against  each 
other.  The  cost  function  for  the  different  initial  values  of  (1  is  shown  in  Figure  3.3. 
It  is  obvious  that  the  algorithm  does  not  converge  to  the  same  solution,  when  it  is 
run  with  different  initial  values. 
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The  initial  value  of  the  spoiler,  ie  ^25  was  varied  between  1.5  and  0.6,  while 
the  initial  value  of  the  elevator  was  held  constant  at  0.5.  The  optimization  history 
is  shown  in  Figure  3.4,  while  Figure  3.5  shows  (2  as  function  of  Ci.  In  this  case  the 
algorithm  converges  to  the  solution  (opt  =  [0.06,0.55]^,  for  all  initial  conditions  but 
the  two  smallest.  Figure  3.6  shows  the  cost  function  for  different  initial  values  of  ^2- 
The  conclusion  is  that  the  calculated  solution  is  sensitive  to  the  initial  value, 
especially  the  initial  value  of  Ci.  The  obvious  conjecture  is  that  the  feasible  set  $  not 
is  convex. 


Figure  3.1:  Optimization  history  of  -  varying  (uniUai 

The  calculated  solution  always  gives  a  very  small  size  of  the  elevator.  That  is, 
however,  consistent  with  the  problem  formulation. 

E.  CONTROLLERS 

The  controllers  for  the  five  flight  conditions  that  are  listed  below,  were  calculated 
for  (  =  [0.1, 0.7]^  which  is  close  to  the  lowest  values  given  in  Figures  (3. 1-3.5).  The 
resulting  five  controllers  are  presented  in  below.  Tables  (3.4-3. 8)  show  the  eigenvalues 


20 
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0.95 


0.9 


0.85 


0.8 


0.75 


Figure  3.2:  Optimization  history:  ^2  function  of  Ci 

and  dampings  for  the  closed  loop  system  for  the  five  flight  conditions.  The  designed 
controllers  give  the  closed  loop  system  rather  slow  dynamics.  This  can  be  corrected 
by  adjusting  the  weights  in  the  C  matrix.  The  damping,  which  varies  between  0.54 
and  0.88  is,  however,  good.  As  the  purpose  of  the  performed  work  is  to  develop  a 
procedure  for  the  design  of  an  controller,  rather  than  actually  designing  a  controller, 
the  results  are  clearly  acceptable. 


■  39.7 

-10.7 

4.2 

29.9 

26.0 

10.0  ■ 

K^  = 

505.3 

-149.3 

0.9 

150.5 

-130.5 

32.7 

-1201.0 

297.1 

-1.7 

-300.7 

258.9 

-65.1  . 

r 

-3.76 

36.86 

17.07 

5.77 

57.94 

2.70 

II 

105.6 

-148.4 

-2.04 

139.9 

-39.57 

31.30 

L 

-154.97 

71.75 

-0.246 

-72.83 

4.70 

-17.23 

15.28 

-2.55 

7.49 

27.44 

40.65 

7.74 

A’®  = 

230.2  - 

-151.4 

0.913 

151.9 

-59.50 

32.15 

-362.5 

153.8  - 

-1.101 

-156.1 

55.94 

-31.67 

38.51 

-19.58 

3.80 

33.03 

29.71 

8.88 

II 

459.4 

-213.9 

1.06 

215.3 

-100.1 

41.96 

-837.7 

295.6 

-1.44 

-298.4 

134.5 

-57.38 
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105.4 

-64.3 

6.8 

87.6 

23.3 

17.5 

787.0 

-505.4 

2.4 

510.5 

-149.9 

89.2 

-1065.3 

543.6 

-2.6 

-550.0 

153.4 

-95.0 

TABLE  3.4:  Closed  loop  eigenvalues  and  damping,  condition  1 


Eigenvalue 

Freq.  (rad/sec) 

-0.2343  -h  0.1751i 

0.8010 

0.2924 

-0.2343  -  0.1751i 

0.2924 

-2.1488 

MMiMiliM 

2.1488 

-2.6423 

1.0000 

2.6423 

-5.2115  -h  4.9987i 

0.7217 

7.2213 

-5.2115  -  4.9987i 

0.7217 

7.2213 

Figures  (3. 7-3. 9)  show  the  actuator  step  response  for  flight  condition  1,  while 
figures  (3.10-3.12)  show  the  Bode  plots  for  the  actuator  dynamics. 
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Figure  3.4:  Optimization  history  of  ^  -  varying  (2initiat 


TABLE  3.5:  Closed  loop  eigenvalues  and  damping,  condition  2 


Eigenvalue 

Damping 

Freq.  (rad/sec) 

-0.2152  +  0.1162i 

0.2446 

0.2446 

-1.4722 

1.4722 

-1.8325  -I-  1.3402i 

-1.8325  -  1.3402i 

-5.7166 

1.0000 

5.7166 

TABLE  3.6:  Closed  loop  eigenvalues  and  damping,  condition  3 


Eigenvalue 

Damping 

Freq.  (rad/sec) 

-0.1905  +  0.1440i 

-0.1905  -  0.1440i 

■ilrfehirM 

-1.3747 

1.0000 

1.3747 

-2.5642 

2.5642 

-3.6290  -1-  3.1582i 

MiMrfegM 

4.8108 

-3.6290  -  3.1582i 

0.7543 

4.8108 
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Figure  3.5:  Optimization  history:  ^2  as  function  of 


0  S  10  IS  20  25  30  35 


Figure  3.6:  Costfunction  for  different  Cimuai 
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Amplitude 


TABLE  3.7:  Closed  loop  eigenvalues  and  damping,  condition  4 


Eigenvalue 


-2.1119 


-3.7694  -K  5.8935i 


-3.7694  -  5.8935i 


Damping  Freq.  (rad/sec) 


0.8 

0.8 


1.0000 


0.5388 


0.5388 
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6.9959 


6.9959 


TABLE  3.8:  Closed  loop  eigenvalues  and  damping,  condition  5 


Eigenvalue 

Damping 

Freq.  (rad/sec) 

-0.1697  -H  0.1127i 

0.8330 

0.2038 

-0.1697  -  0.1127i 

0.8330 

0.2038 

-2.0966 

1.0000 

2.0966 

-3.4113 

1. 

0000 

3.4113 

-4.3355  -H  3.9329i 

7407 

5.8536 

-4.3355  -  3.9329i 

7407 

5.8536 

Input  1  Output  1 


Input  1  Output  2 


Tirne  (socg) 
Input  1  Output  3 


Time  (eece) 
Input  2  Output  1 


Time  (secs) 


Time  (secs) 


Figure  3.7:  Actuator  step  response 


Frequency  (rad/sec) 


Figure  3.11:  Bode  plot  for  spoiler 


IV.  CONCLUSIONS  AND 
RECOMMENDATIONS 


A.  CONCLUSIONS 

This  thesis  develops  the  technique  of  integrating  the  design  of  aircraft  plant 
parameters  and  the  design  of  the  corresponding  feedback  controller.  The  developed 
algorithm  using  Matlab  and  LMTLab  gives  a  solution  to  the  problem  of  minimizing 
the  control  surfaces  of  an  aircraft,  while  finding  a  feasible  controller.  The  result, 
when  this  algorithm  is  applied  to  the  longitudinal  dynamics  of  an  F-4  aircraft,  is  a 
set  of  controllers,  one  for  each  flight  condition,  and  a  single  value  of  the  sizes  of  the 
moving  horizontal  tail  and  the  spoiler. 

The  approach  used  to  formulate  an  integrated  plant  optimization  and  Hoo  con¬ 
troller  design  as  a  linear  matrix  inequality  results  in  a  non-convex  optimization  prob¬ 
lem. 

Compared  to  [Ref.  1]  the  problem  formulation  has  been  expanded  to  handle 
multiple  flight  conditions  and  to  include  thrust  control.  As  the  developed  algorithm 
using  LMI-Lab  is  more  computationally  efflcient  than  the  previous  algorithm,  it  is 
possible  to  handle  much  larger  problems  than  in  [Ref.  1]. 

The  calculated  solutions  do,  however,  depend  on  the  initial  value,  as  the  problem 
is  non-convex.  It  is  thus  necessary  to  run  the  algorithm  several  times,  to  be  reasonably 
sure  that  a  good  quality  solution  has  been  obtained. 

This  work  has  increased  the  knowledge  on  how  integrated  plant  optimization 
and  controller  design  problems  should  be  formulated  and  provides  a  good  ground  for 
further  work.  The  completion  of  the  recommended  future  work  in  the  next  section  will 
maJce  this  approach  to  integrated  plat/controller  design  a  useful  tool  for  engineers. 
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B.  RECOMMENDATIONS 


The  convergence  properties  of  the  optimization  algorithm  are  unclear  and  should 
be  investigated.  The  first  issue  to  consider  is  whether  the  algorithm  converges  to  the 
global  minimum  for  a  convex  problem.  The  second  issue  is  to  investigate  the  possi¬ 
bilities  of  finding  a  convex  hull  to  the  feasible  set  of  plant  parameters  and  controllers, 
and  establish  a  lower  bound  for  the  optimal  solution. 

The  problem  treated  in  this  thesis  only  concerns  the  dynamics  of  the  aircraft 
for  small  perturbations  around  a  trimmed  equilibrium.  The  formulation  does  not 
include  the  fact  that  the  aircraft  configuration  given  by  the  optimal  value  of  (  must 
be  possible  to  trim  at  the  given  flight  condition.  The  problem  formulation  must  be 
changed  to  include  the  requirement  that  the  aircrft  must  be  possible  to  trim. 

The  performance  requirements  concerning  maneuverability  are  translated  to  an 
Hoo  formulation  based  on  the  linearized  model  developed  for  small  perturbations 
around  an  equilibrium.  A  non-linear  dynamical  model  should  instead  be  used  for  the 
maneurability  requirements. 

The  result  of  this  work  is  a  set  of  time  invariant  feedback  controllers,  each 
of  them  valid  for  one  flight  condition.  This  set  of  controllers  can  be  used  in  the 
construction  of  a  gain  scheduled  controller,  that  works  for  the  non-linear  aircraft 
dynamics  at  all  flight  conditions.  The  resulting  closed  loop  system  is  time  variant 
and  non-linear.  The  gain  scheduling  must  be  done  such  that  when  this  system  is 
linearized  at  different  flight  conditions,  the  linear  system  has  the  desired  properties. 
The  linearized  system  dynamics  should  be  as  independent  of  the  gain  scheduling 
parameter  as  possible.  A  continuation  of  the  work  presented  in  this  thesis  is  to  find 
a  way  to  interpolate  between  the  different  controller  such  that  this  is  achieved. 


30 


APPENDIX  A:  MATLAB  CODE 


In  the  appendix  the  Matlab  code  using  LMI-Lab  is  listed.  The  main  program  is 
f4opt.m,  which  uses  the  two  functions  sol4.m  and  solS.m.  This  combination  of  Matlab 
and  LMTLab  needs  to  be  run  on  a  SUN  Sparc  station  and  does  not  work  on  other 
computers. 

In  f4opt.m  the  LMI’s  are  formulated  using  two  equivalent  versions  of  the  LMI 

■  A‘(C)T‘  +  +  BiiQW'  +  (C)  B{{C)  (CV*  +  DW'f 

R{W\Y\0=  Bf(C)  -I  0 

C'F’  +  D^W^  0  -I 

these  are 

■  A*(C)T*'  +  +  BiiC)W^  +  (0+ 

i21(IT‘,y‘,0=  +51(05^(0  {C^Y'  +  D^W')^  < 

C'Y^  +  DW'  -I 


that  is  used  for  solving  for  the  controller,  i.e.  W  and  Y,  and 


R2{W\Y\C) 


yi'(c)y*  +  yM‘^(C)  +  B'^iow'  +  (0+ 

+{C'Y^  +  D^W'fiC^Y^  + 

BUY 


that  is  used  to  solve  for  the  optimal  C. 
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•/.  This  program  performs  an  integrated  minimization  of  elevator  and  spoiler, 

•/,  and  calculation  of  H-infinity  controller.  The  aircraft  is  an  F-4  with 
%  aircraft  data  from  Schmidt.  The  calculations  are  performed  for  five 
*/.  different  flight  conditions. 

^  The  program  uses  the  toolbox  to  be  LMI~Lab  and  the  two  matlab  functions 

y,  sol4  and  sol5. 

y,  The  system  is  modeled  as 

y,  x_dot=A*x+B_l*w+B_2*u 

y,  z=C*x+D_2*u 

y,  The  program  calculates  the  H-infinity  controller  K=W*inv(Y)  (u=K*x) 

y,  To  use  this  program  to  design  a  controller  the  weights  in  the  B_l,  B_2,  C  and 

y,  D_2  matrices  can  be  changed. 

y,  The  output  of  the  program  is  stored  in  the  arrays  lam  <maximum  eigen  value 
y,  for  each  LMI),  zet  (plant  parameters)  and  cost  (cost  function). 

clear; 

g*32.174;  y.  ft/s“2 
m=1210;  y,  slugs 
c=16.0;  y.  ft 
S=530;  y,  ft “2 

Xdlc=-0.02*g*57.1/40;  '/,  drag  from  spoiler 
Zdlc=0.2*g*57. 1/40;  '/,  destroyed  lift  from  spoiler 
c=[l;l];  y,  weights  in  cost  function 

Flight  condition  6  ****** 

Iy=122190;  X  slug*ft 
U*1452;  y.  ft/s 

Xu=-0.0072;  '/,  dimensional  aerodynamic  derivatives 
Xa*-29.98; 

Zu*-0.515; 

Za=-716.7; 

Zad*-0.519; 

Zq=-2.24; 

Mu=0.0025; 

Ma=-28.94; 

Had»-0.122; 

Mq=-0.488; 

Xd=0.00; 

Zd=-70.67; 

Md*-16.00; 

xc=Md/Zd*Iy/m;  */,  lever  for  elevator 
alpha0=2.6/57.3;  */,  angle  of  attack  in  radians 
theta0*alpha0 ; 
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thrust _acc=0.05*g;  '/,  max  acceleration  from  thrust 


sigma_g®l/57.3;  */,  max  commanded  flight  path  angle 

sigma. a=5/U;  '/,  wind  gusts 

Umax  1=20/57. 3;  */,  max  elevator  deflection 

umax2=40/57.3;  */,  max  spoiler  deflection 

sigma. aout=l .5/57.3;  %  max  error  in  alpha 

C=[zeros(3,6) ;0  1/sigma.aout  000  0;0  0  0  0  60  0;0  0000  10]; 
D2=[l/umaxl  0  0;0  l/umax2  0;0  0  1/thrust. acc;zeros(3, 3)] ; 

*/,  A  matrix  for  nominal  plant 

An=[U*Xu  Xa  0  -g*cos(theta0) ;  U*Zu  Za  U+Zq  -g*sin(theta0) ; 

U*Mu  Ma  Mq  0;  0010]; 

In=[U  0  0  0;0  U-Zad  0  0;0  -Mad  10;  0001]; 

A006=[inv(In)*An  zeros(4,2);0  1  0  0  0  0;0  -1  0  1  0  0] ; 

'/,  Calculation  of  derivatives  for  wing/body 
Za=Za~Zd; 

Ma=Ma-xc*Zd*m/Iy ; 

Zq=Zq-xc*Zd/U ; 

Mq=Hq-xc*xc*Zd*m/Iy/U ; 

Zad=Zad-xc*Zd/U ; 

Mad=Mad-xc*xc*Zd*m/Iy/U ; 

An=[U*Xu  Xa  0  -g*cos(theta0) ;  U*Zu  Za  U+Zq  -g*sin(theta0) ; 

Ut'Mu  Ma  Mq  0;  0  0  10]; 

invIn=[l/U  0  0  0;0  1/(U-Zad)  0  0;0  Mad/(U-Zad)  1  0;  0  0  0  1]; 
A06=[invIn*An  zeros(4,2);0  1000  0;0  -10100]; 

dlnd=[0  0  0  0;0  1  0  0;  0  Mad+xc*xc*m/Iy*(U-Zad)  0  0;0  0  0  0]/U/(U-Zad)“2; 
dAnd=[0  0  0  0;0  1  1/U  0;  0  xc*m/Iy  xc*xc*m/Iy/U  0;0  000]; 

A16=[(dInd*An+inv(In)*dAnd)*Zd  zeros(4,2);  zeros(2,6)]; 


B106=zeros(6,2);  B116=B106;  B126=B106; 

B106(l;4,l)=A06(l:4.2)*sigma.a; 

B106(6 , 2)=sigma.g; 

dB106d=[0  0;1  0;xc*m/ly  0;0  0]*sigma_a; 
B116=[(dInd*B106(l:4, :)+invIn*dB106d)*Zd;zeros(2,2)]  ; 

B206=[0  0  l/U;zeros(5,3)] ; 

dB206d=[0  0  0;1  0  0;xc*m/ly  0  0;0  0  0]; 
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dB206s=[0  Xdlc/Zdlc  0;0  1  0;0  0  0;0  0  0]; 

B216=[(dInd*B206(l :4, : )+invIn*dB206d)*Zd;  zeros(2,3)] ; 
B226=[invIn*dB206s*Zdlc;  zeros (2, 3)] ; 

Flight  COndltlOIl  5  **♦*♦***♦*******♦*♦**♦*♦* 

U=1167; 

Xu=-0.0136; 

Xa=-16.53; 

Zu=-0.747; 

Za=-848.3; 

Zad=-1.24; 

Zq=-4.09; 

Mu=0.0022; 

Ma=-29.03; 

Mad»-0.288; 

Mq=-0.746; 

Xd=0.00; 

Zd=-90.44; 

Md=-20,70; 
xc*Md/Zd*Iy/m; 
alpha0=1.6/57.3; 
thGtaO=alphaO ; 
sigma_a*5/U; 

An=[U*Xu  Xa  0  -g*cos(thetaO) ;  UkZu  Za  U+Zq  -g*sin(thetaO) ; 

U*Mu  Ha  Mq  0;  0  0  1  0] ; 

iiivIn»Cl/U  0  0  0;0  1/(U-Zad)  0  0;0  Mad/(U-Zad)  10;  0  0  0  1]; 

A05= [invIn*An  zeros(4,2);0  1000  0;0  -10100]; 

dlnd=[0  000;0100;0  Mad+xc*xc*m/Iy*(U-Zad)  0  0;0  0  0  0]/U/(U-Zad)“2; 
dAnd=[0  0  0  0;0  1  1/U  0;  0  xc*m/Iy  xc*xc*in/Iy/U  0;0  0  0  0]; 

A15=[(dInd*An+invIn*dAnd)*Zd  zGros(4,2);  zeros(2,6)]; 

B105*zeros(6,2);  B115*B105;  B125*B105; 

B105(l:4,l)=A05(l:4,2)*sigma.a; 

B105(6,2)=sigma_g; 

dB105d»C0  0;1  0;xc*m/ly  0;0  0]*sigma_a; 

B115=[(dInd*B105(l:4, :)+invIn*dB105d)*Zd;zeros(2,2)] ; 

B205=[0  0  l/U;zeros(5,3)] ; 

dB205d=[0  0  0;1  0  0;xc*m/ly  0  0;0  0  0]; 

dB205s=C0  Xdlc/Zdlc  0;0  1  0;0  0  0;0  0  0]; 

B215=[(dInd*B205(l:4, :)+invIn*dB205d)>fZd;  zeros(2,3)] ; 
B225=[invIn*dB20Ss*Zdlc;  zeros (2, 3)] ; 
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Flight  condition  4  ***************************** 


U=867; 

Xu=-0.0123: 

Xa=-7.43; 

Zu=-0.571; 

Za=-475.4; 

Zad=-1.01; 

Zq=-2.89; 

Mu=-0.0028; 

Ma=-7.877; 

Mad= -0.234; 

Mq=-0.487; 

Xd=0.00: 

Zd=-49.65; 

Md=-11.40; 

xc=Md/Zd*Iy/m; 

alpha0=2.6/57.3; 

thetaO=alphaO; 

sigma_a=5/U; 


An®[U*Xu  Xa  0  -g*cos(thetaO) ;  U*Zu  Za  U+Zq  -g*sin(thetaO) ; 

U*Mu  Ha  Hq  0;  0  0  1  0] ; 

invIn=[l/U  0  0  0;0  1/(U-Zad)  0  0;0  Mad/(U-Zad)  1  0;  0  0  0  1]; 

A04= [invIn*An  zeros(4,2);0  1000  0;0  -10100]; 

dlnd=[0  000;0100;0  Mad+xc*xc*m/Iy*(U-Zad)  0  0;0  0  0  0]/U/(U-Zad)‘2; 
dAnd=C0  0  0  0;0  1  1/U  0;  0  xc*m/Iy  xc*xc*m/Iy/U  0;0  0  0  0]; 

A14=[(dInd*An+invIn>t<dAnd)*Zd  zeros(4,2);  zeros(2,6)]; 

B104=zeros(6,2);  8114=8104;  8124=8104; 

8104(1 :4,l)=A04(l:4,2)t<sigma_a; 

8104(6,2)=sigma_g; 

dB104d=[0  0;1  0;xc*m/ly  0;0  0]*sigma_a; 

B114=[(dInd*B104(l :4, : )+invIn*dB104d)*Zd;zeros(2,2)]  ; 

B204=[0  0  l/U;zeros(5,3)] ; 

dB204d=[0  0  0;1  0  0;xc*m/ly  0  0;0  0  0]; 

dB2048=C0  Xdlc/Zdlc  0;0  1  0;0  0  0;0  0  0]; 

B214=[(dInd*B204(l :4, : )+invIn*dB204d)*Zd;  zeros(2,3)] ; 

8224= [invIn*dB204s*Zdlc ;  zeros (2,3)] ; 

Flight  COnditiOIl  3 
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U=584; 

Xu*-0.0176; 

Xa=-24.25; 

Zu=-0.337; 

Za=-162.2; 

Zad=-0.591; 

Zq=-1.82; 

Mu=-0.0001; 

Ma=- 1.927; 

Had=-0.141; 

Mq=-0.307; 

Xd=0.00; 

Zd=-20.98; 

Hd=-4.90; 

xc=Md/Zd*Iy/m; 

alpha0=9.4/57.3; 

thetaO=alphaO; 

sigma_a=5/U; 

An=[U*Xu  Xa  0  -g*cos(thetaO) ;  U*Zu  Za  U+Zq  -g*siii(thetaO) ; 

U*Mu  Ma  Mq  0;  0  0  1  0] ; 

invIn=[l/U  0  0  0;0  1/(U-Zad)  0  0;0  Mad/(U-Zad)  10;  0001]; 

A03= [invIn*An  zeros(4,2);0  1000  0;0  -10100]; 

dlnd=[0  000;0100;0  Mad+xc*xc*m/Iy*(U-Zad)  0  0;0  0  0  0]/U/(U-Zad)"2; 
dAnd=[0  0  0  0;0  1  1/U  0;  0  xc*m/Iy  xc*xc*m/Iy/U  0;0  0  0  0]; 

A13=[(dInd>t'An+invIn*dAnd)*Zd  zeros(4,2);  zeros(2,6)]; 

B103=zeros(6,2) ;  B113*B103;  B123=B103; 

B103(1:4,1)=A04(1 :4,2)*sigma_a; 

B103(6,2)=sigma_g; 

dB103d=[0  0;1  0;xc*m/ly  0;0  0]*sigma_a; 

B113=C(dInd*B103(l:4, :)+invIn*dB103d)*Zd;zeros(2,2)] ; 

B203=[0  0  l/U;zeros(5,3)] ; 

dB203d=C0  0  0;1  0  0;xc*m/ly  0  0;0  0  0]; 

dB203s=C0  Xdlc/Zdlc  0;0  1  0;0  0  0;0  0  0]; 

B213=[(dInd*B203(l :4, :)+invIn*dB203d)*Zd;  zeros(2,3)] ; 

B223= [invIn*dB203s*Zdlc ;  zeros(2,3)] ; 


Flight  condition  2  *********************** 
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U=952; 

Xu=-0.0215; 

Xa=-5.93; 

Zu=- 1.170; 

Za=-1103.2; 

Zad=-2.00; 

Zq=-6.00; 

Hu=-0.0044: 

Ma*-17.00; 

Mad=-0.457; 

Mq=-0.993; 

Xd=0.00; 

Zd=-107.0; 

Md=-25.00; 
xc=Md/Zd*Iy/in; 
alpha0*0.5/57.3; 
theta0=alpha0 ; 
sigma_a=5/U; 

An=[U*Xu  Xa  0  -g+cosCthetaO) ;  U*Zu  Za  U+Zq  -g*sin(thetaO) ; 

U*Mu  Ma  Mq  0;  0  0  1  0] ; 

invIn=[l/U  0  0  0;0  1/(U-Zad)  0  0;0  Had/(U-Zad)  10;  0001]; 

A02= [invIn*An  zeros(4,2);0  1000  0;0  -10100]; 

dlnd=[0  000;0100;0  Mad+xc*xc*m/Iy*(U-Zad)  0  0;0  0  0  0] /U/(U-Zad)“2 
dAnd=C0  0  0  0;0  1  1/U  0;  0  xc*m/Iy  xc*xc*m/Iy/U  0;0  000]; 

A12=C(dInd*An+invIn*dAnd)*Zd  zeros(4,2);  zeros(2,6)]; 

B102=zeros(6,2);  B112=B102;  B122=B102; 

B102(l:4,l)=A02(l:4,2)*sigma_a; 

B102(6,2)=sigma_g; 

dB102d*[0  0;1  0;xc*m/ly  0;0  0]*sigma_a; 
B112=C(dInd*B102(l:4,:)+invIn*dB102d)*Zd;zeros(2,2)] ; 

B202=C0  0  l/U;zeros(&.3)] ; 

dB202d*[0  0  0;1  0  0;xc*m/ly  0  0;0  0  0]; 

dB202s=C0  Xdlc/Zdlc  0;0  1  0;0  0  0;0  0  0]; 

B212=C(dInd*B202(l:4,:)+invIn*dB202d)*Zd;  zeros (2, 3)] ; 

B222*CinvIn*dB202s*Zdlc;  zeros (2, 3)] ; 

^  )(C  *  it  *  *  it  }|(  )((***********  )|(  *  >|(  *  t  )|(  ********  ^  *)*(*** 

y,  Calculations 

^i(ieici(ici(ici(ici(iticic******************  ******************************** 

'/.  Inital  value  of  zeta 
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cost=  []  ; 
zet-  []  ; 


zeta*[0. 1;0.7] ; 
costold®c*zeta ; 

for  iter® 1:40 

A2®A02+A12*zeta(l) ; 

B12=B102+B112*zeta(l) ; 

B22=B202+B212*zeta(l)+B222*zeta(2) ; 

A3=A03+A13*zeta(l) ; 

B13=B103+B113*zeta(l) ; 

B23=B203+B213*zeta(l)+B223*zeta(2) ; 

A4=A04+A14*zeta(l) ; 

B14=B104+B114*zeta(l); 

B24=B204+B214*ZGta(l)+B224*zeta(2) ; 

A5=A05+A15*zeta(l): 

B15=B105+B115*zeta(l) ; 

B25®B205+B215*zeta(l)+B225*zeta(2); 

A6=A06+A16*zeta(l): 

B16=B106+B116*zeta(l) ; 

B26=B206+B216>»zeta(l)+B226*zeta(2) ; 

*/,  Solve  LMI’s  for  controller  in  iteration  1  sol  4  is  used,  whil  sol5  is 
y,  used  in  the  following  iterations  to  correct  the  solution  when  feasp 
y,  gives  a  non-negative  definite  solution.  In  that  case  the  old  solution 
y,  is  retained. 

if  iter<2 

[Y2  H2  maxe2  x2]=sol4(A2,B12,B22,C,D2) ; 

[Y3  W3  maxes  x3]=sol4(A3,B13,B23,C,D2) ; 

[Y4  W4  maxe4  x4]=sol4(A4,B14.B24,C,D2) ; 

[Y5  W5  maxes  x5]=sol4(A5,B15,B25,C,D2) ; 

[Y6  W6  maxes  x6]=sol4(A6,B16,B26,C,D2) ; 
else 

[Y2  W2  maxe2  x2]®sol5(A2,B12,B22,C,D2,x2) ; 

[Y3  W3  maxes  x3]=sol5(A3,B13,B23,C,D2,x3) ; 

[Y4  W4  maxe4  x4]=sol5(A4,B14,B24,C,D2,x4) ; 

[Y5  W5  maxes  x5]=solS(A5,B15,B2S,C,D2,x5) ; 

[Y6  W6  maxes  x6]=sol5(A6,B16,B26,C,D2,x6) ; 
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end; 

ine=[inaxe2  maxeS  inaxe4  maxeS  maxe6]  ; 
lam=[lam;  me]  ; 

'/.Solve  LMI’s  for  optimal  zeta 

R112=A02*Y2+Y2*A02 ’ +B202*W2+W2 ’ *B202 ’ + (C*Y2+D2*W2) ’ * (C*Y2+D2*W2) 
R1 13=A03*Y3+Y3*A03 ’ +B203*W3+W3 ’ *B203 ’ + (C*Y3+D2*W3) ’ * (C*Y3+D2*W3) 
R1 14=A04*Y4+Y4*A04 ’ +B204*W4+W4 ’ *B204 ’ + (C*Y4+D2*W4) ’ * (C*Y4+D2*W4) 
R1 15=A05*Y5+Y5*A05 ’ +B205*W5+W5 ’ *B205 ’ + (C*Y5+D2*W5) ’ * (C*Y5+D2*W5) 
R1 16=A06*Y6+Y6*A06 ’ +B206*W6+W6 ’ *B206 ’ + (C*Y6+D2*W6) ’ * (C*Y6+D2*W6) 


lm=  []  ; 

lm=add_lmi(lm, [6  2] ,8) 
lm®add_lmi(lm, [6  2] ,8) 
lm=add_lmi(lm, [6  2]  ,8) 
lm*add_lmi(lm, [6  2]  ,8) 
lm=add_ Imi ( Im , [6  2] , 8 ) 

lm=add_lmi(lm, [1] ,1) ; 
lm=add_lmi(lm, [1] , 1) ; 


lm=add_var(lm,  1,  [1  1]); 
lm=add_var(lm,  1,  [11]); 


lm*add_term(lm, 

[1 

1 

1 

0], 

R112); 

lm=add_t  erm ( Im , 

[1 

1 

1 

1], 

B212.  W2, 

’s’) ; 

lm=add_term(lm. 

[1 

1 

1 

1], 

A12.  Y2, 

’s’); 

lm=add_term(lm, 

[1 

1 

1 

2], 

B222,  W2, 

’s’); 

lm=add_term(lm. 

[1 

1 

2 

0], 

B102) ; 

lm=add_term(lm. 

[1 

1 

2 

1], 

B112) ; 

lm=add_term(lm. 

[1 

2 

2 

0], 

-1): 

lm=add_term(lm. 

[2 

1 

1 

0], 

R113); 

lm=add_term(lm, 

[2 

1 

1 

1], 

B213,  W3, 

’s’); 

lm=add_term(lm. 

[2 

1 

1 

1], 

A13,  Y3, 

’s’); 

lm®add_term(lm, 

[2 

1 

1 

2], 

B223,  W3, 

’s’); 

lm=add_term(lm , 

[2 

1 

2 

0], 

B103); 

lm®add_t erm ( Im , 

[2 

1 

2 

1], 

B113); 

lm*add_term(lm, 

[2 

2 

2 

0], 

-1): 

lm*add_term(lm, 

[3 

1 

1 

0], 

R114) ; 

lm=add_t  erm ( Im , 

[3 

1 

1 

1]. 

B214,  W4, 

’s’); 

lm=add_term(lm, 

[3 

1 

1 

1], 

A14,  Y4. 

’s’); 

lm=add_term(lm. 

[3 

1 

1 

2], 

B224,  W4, 

’s’); 
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lm=add_term(lin, 

[3 

1 

2 

03, 

B104) ; 

lm=add_term(lm, 

[3 

1 

2 

13. 

B114); 

lin=add_term(lm , 

[3 

2 

2 

03. 

-1); 

lm=add_term(lm , 

[4 

1 

1 

03. 

R115): 

lm=add_term(lm, 

C4 

1 

1 

1]. 

B215,  W5, 

lin=add_term  (Im , 

[4 

1 

1 

13. 

A15,  Y5, 

lm=add_t erm ( Im , 

[4 

1 

1 

23. 

B225,  W5. 

lm=add_tenn(lm, 

[4 

1 

2 

03. 

B105); 

lm=add_term (Im , 

[4 

1 

2 

13. 

B115); 

lm=add_term(lm , 

[4 

2 

2 

03. 

-1); 

lm=add_term(lm, 

[5 

1 

1 

03. 

R116); 

lm=add_term(lm, 

[5 

1 

1 

13. 

B216,  W6. 

lm=add.term(lm, 

[5 

1 

1 

13. 

A16,  Y6, 

lm=add_term(lm, 

[5 

1 

1 

23, 

B226,  W6, 

lm=add_term(lm, 

[5 

1 

2 

03. 

B106); 

lm=add_term(lm, 

[6 

1 

2 

13. 

B116); 

lm=add_term(lm, 

[5 

2 

2 

03. 

-1); 

lm=add_t erm (Im , 

[-e 

;  1 

L  j 

L  13. 

1.  1); 

lm=add_term(lm, 

[-7 

'  j 

L  1 

L  23. 

1,  1); 

Ccopt,xopt3=linobj (lm,c, [0  0  0  50  1] ,zeta) ; 

iter 

zeta 

^Controllers  for  each  flight  condition 
K2n=W2*inv(Y2); 

K3n=W3*inv(Y3); 

K4n=W4*inv(Y4) ; 

K5n=WS*inv(Y5); 

K6n=W6*inv(Y6) ; 

K2=[K2;  K2n]; 

K3=CK3;  K3n]; 

K4=[K4;  K4n3; 

K5=CK5;  K5n3; 

K6=CK6;  K6n3; 

Cost=c*xopt 
cost=[cost;  Cost3 ; 
zet=[zet;  zeta’ 3; 
xopt 
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zeta=xopt ; 

if  abs (Cost-Costold) / Costold<0 . 0001 

if  max (me) <0 

break; 

end; 

end; 

'/,  Closed  loop  system  in  the  final  version 
Ac2=A2+B22*K2n; 

Ac3*A3+B23*K3n; 

Ac4=A4+B24*K4n; 

Ac5=A5+B25*K5n; 

Ac6=A6+B26*K6n ; 

damp(Ac2} 

damp(Ac3) 

damp(Ac4) 

damp(Ac5) 

damp(Ac6) 

D=zeros(3,3) ; 

end 
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function  [Y,  W,  maxe,  x]=sol4(A,Bl,B2,C,D) 
y.  Solves  the  H-inf  control  problem  and  returns  Y,  W  matrices 
y,  the  maximum  eigenvalue  euid  the  solution  vector  in  LMI  format 
Im®  [] ; 

lm=add_lmi(lm, [6  6] , 12) ; 
lm=add_lmi(lm, [6] ,6) ; 
lm=add_var(lm,  1,  [61]); 
lm=add_var(lm,  2,  [3  6]); 
lm=add_term(lm, [1  110],  B1*B1’); 
lm=add_term(lm,  [1  1  1  1],  A,  eye (6),  ’s’); 
lm=add_term(lm,  [1  1  1  2],  B2,  eye(6),  ’s’); 
lm®add_term(lm,  [1  2  1  1] ,  C) ; 
lm=add_term(lm,  [1  2  1  2],  D); 
lm=add_term(lm,  [1220],  -1); 
lm=add_term(lm,  [-2  111],  eye(6),  eye(6)); 

[tmin,xf eas]®feasp(lm, [200  0  100  1],0); 
evalsys=eval_lmi(lm, [xfeas] ) ; 

[lhs,rhs]=show_lmi(evalsys, 1) ; 

Y»dec2mat(lm,xfeas,l) ; 

W®dec2mat(lm,xf eas,2) ; 
maxe=m€oc  (eig(lhs)) 
x=xfeas; 
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function  [Y,  W,  maxe,  x]=sol5(A,Bl»B2,C,D,xl) 

y,  Solves  the  H-inf  control  problem  given  A,  Bl,  B2,  C,  D  matrices 

y,  and  the  previous  LMI  solution.  This  function  checks  if  all  the 

'/•  eigenvalues  are  negative  and  if  they  are  not  returns  a  solution  based 

y,  on  the  previous  iteration. 

lm=  []  ; 


Im* 

lm= 

lm= 

lm= 

lm= 

Im* 

lm= 

Im* 

lm= 

lm= 

lm= 


add_lmi(lm,[6  6], 12); 
add_lmi(lm, [6] ,6) ; 
add_var(lm,  1,  [61]); 
add_var(lm,  2,  [36]); 
add.termClm, [1  110],  B1*B1’); 
add_term(lm,  [1111] 

[1112] 


add_term(lm, 

add_term(lm, 

add.termClm, 

add.termClm, 

add.termClm, 


[1 

[1 

[1 

[-2 


1 
1 
1 
2 

1  1  1] 


1] 

2] 

0] 


A,  eyeC6),  ’s’); 
B2,  eyeC6),  ’s’); 

C) ; 

D) ; 

-1); 

eyeC6),  eyeC6)); 


6valsys*eval.lmiClm, [xl] ) ; 
[lhs,rhs]=show.lmiCevalsys,l) ; 
maxe 1 =max C e ig C Ihs ) ) 

[tmin,xfeas]=f easpClm, [200  0  100  1],0); 
evalsys=eval.lmiClm, [xf eas] ) ; 

[Ihs ,rhs]=show.lmiCGvalsys , 1) ; 
maxeamaxCeigClhs) ) 


if  maxe<0 
x=xf eas ; 
else 

maxe=maxel ; 
x=xl ; 
end; 


Y=dec2matClni,x,l) ; 
W=dec2matClm,x,2) ; 
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